DAMPING OF DROP OSCILLATIONS BY SURFACTANTS AND SURFACE VISCOSITY 

Brian M. Rush, Ali Nadim, Dept. Aero, and Mech. Engr., Boston University, Boston, MA 02215 USA 


Abstract 

The roles played by surfactants and interfacial rheology 
on damping the shape oscillations of liquid drops are 
analyzed for the case of axi symmetric shape oscillations 
of a nearly spherical liquid drop surrounded by another 
fluid in the absence of gravity. Both fluids are taken to 
be viscous, although the Reynolds number associated 
with the shape oscillations is assumed large enough that 
deviations from inviscid flow are only present in thin 
Stokes boundary layers near the no-slip interface. Also, 
an insoluble surfactant is assumed to be present at the 
interface and surface tension is taken to be a linearly de- 
creasing function of local surfactant concentration. This 
surfactant layer is further assumed to behave as a two- 
dimensional Newtonian fluid layer characterized by sur- 
face shear and dilatational viscosities. 

Under these conditions, several sources can be iden- 
tified for mechanical dissipation and the ultimate damp- 
ing of the shape oscillations of the drop. These include 
viscous effects associated with the bulk fluids that ap- 
pear in two distinct forms: one associated with the os- 
cillatory viscous boundary layers which form near the 
interface between the drop and surrounding fluid, and 
the other associated with the flows far from the inter- 
face which resemble potential flow although they dissi- 
pate energy through weak viscous action. Surfactants 
and surface rheology additionally contribute to dissipa- 
tion in other ways. The surface shear and dilatational 
viscosities dissipate energy within the interface in much 
the same way as the viscous dissipation in the three- 
dimensional bulk fluids just mentioned. Moreover, as 
various parts of the interface locally increase or decrease 
in area during shape oscillations, the concentration of 
surfactants locally decreases or increases. This leads 
to gradients in surfactant concentration on the interface 
where the process of gradient diffusion of surfactants 
within the interface, itself an irreversible process, leads 
to energy dissipation. Also, the Marangoni flow re- 
sulting from this non-uniformity in surface tension con- 
tributes to viscous damping. 

This paper outlines the derivation of a general me- 
chanical energy equation for such a system. It con- 
tains dissipation terms accounting for each of the mech- 
anisms described above. The energy equation is applied 
to a slightly perturbed axisymmetric drop oscillating in 
a low-density surrounding fluid to derive an approxi- 


mate ordinary differential equation (resembling that of 
a damped harmonic oscillator) which describes the time 
evolution of pure shape modes. 

In parallel to the analytical treatment, the imple- 
mentation of a boundary integral method for potential 
(i.e. inviscid) flow is presented for the case of a two- 
dimensional drop oscillating in vacuum. The effect of a 
constant surface dilatational viscosity is included in the 
computations by combining the tangential and normal 
components of the dynamic boundary condition into a 
single equivalent expression. This expression, combined 
with Bernoulli's equation for the pressure, the kinematic 
boundary condition and the regularized Fredholm inte- 
gral equation of the second kind representing Laplace's 
equation for potential flow, produces a coupled set of 
nonlinear equations that allow the time evolution of the 
drop's shape to be followed. Surface dilatational vis- 
cosity is shown to have a damping effect on the free 
oscillations of the drop. 


1 Introduction 

Shape oscillations of drops and bubbles have been stud- 
ied since the works of Kelvin [6] and Rayleigh [15] who 
determined the linearized frequencies of inviscid shape 
modes. Their work was extended by Lamb [7], Reid 
[16] and Valentine et al [19], who included estimates 
of damping by weak viscous effects in the bulk. Miller 
& Scriven [1 1] and Marston [10] further identified the 
important role played by the boundary layers near the 
interface in damping the oscillations. The additional ef- 
fects of surfactants and surface rheology, known to have 
a strong influence on the base frequency and damping 
rates of drops [1] and bubbles [3], have been analyzed 
by Lu & Apfel [8] and Tian et al [18]. Numerical 
studies, based upon the boundary integral method, of 
the dynamics of weakly viscous drops were initiated by 
Lundgren and Mansour [9] and have been extended to 
include surfactant effects by Tian et al Apfel et al [1] 
used numerical simulations in conjunction with experi- 
mental studies of drops in microgravity to quantify the 
important role of surfactants in such systems. 

The present contribution also focuses on the role of 
surfactants and surface rheology on drop oscillations in 
the absence of gravity. In addition to Marangoni effects 
which arise due to non-uniformity of surface tension 


615 



2 THE ENERGY EQUATION 


along the interface in the presence of surfactants, the in- 
terfacial layer of adsorbed surfactants may also possess 
separate rheological properties [4] which change the dy- 
namics. Our goal is to identify the specific mechanisms 
through which surfactants and surface rheology affect 
the system and quantify each of them. An approach 
based on a global mechanical energy balance is outlined 
and allows the damping rates of pure shape modes by 
bulk viscosity, surface viscosity, boundary layer dissi- 
pation, and surfactant transport to be quantified. This 
energy equation generalizes the work of Hsu & Apfel 
[5], who used a simplified energy equation to approxi- 
mate the rate of damping of a viscous drop with a con- 
stant surface tension oscillating in the quadrupole mode. 
Supplementing the analytical treatment, a numerical im- 
plementation of the boundary integral method for poten- 
tial flow is described which incorporates the effects of 
surface viscosity. The method is presented for the os- 
cillations of two-dimensional liquid drops possessing a 
constant surface dilatational viscosity. 


2 The Energy Equation 

Consider a liquid drop of density p and viscosity p to 
be suspended in an infinite medium of density p and 
viscosity //, in the absence of gravity. If both fluids arc 
assumed to be incompressible and Newtonian, the con- 
tinuity and momentum equations take the forms: 

By 

v • v = o , p— = v n for xer(f), (i) 

Dx - - 

V-v = 0, ^ = for (2) 

where v and v respectively refer to the medium and 
drop velocity fields, Y(t) and V(t) are the material vol- 
umes of the medium and drop, and the stress tensors II 
and ft are given by 

II = —pi + 2/tE , E= I[(Vv) + (Vv) r ],(3) 

fl = -pi + 2/7E , E = I[(Vv) + (Vv) 7 ]. (4) 

Here p and p represent the pressures in the two fluids, 
I is the isotropic unit tensor, and E and E are the sym- 
metric and traceless rate-of-strain tensors. 

These field equations need to be supplemented by 
boundary conditions at infinity - that the velocity field 
vanishes and the pressure tends to a constant value - 
and at the interface 5(f) between the drop and medium. 
The interface is assumed to be covered with surfactants 


and therefore possesses its own rheological properties, 
which may be characterized by the surface stress tensor 
U s [4]. The boundary conditions at the interface thus 
assume the respective forms [4]: 

= V * 1 for x e 5(f). 

n (n - n) - -v, n* / w 

(5) 

The velocity at the interface is denoted by v s and is 
equal to the fluid velocities in the medium and drop 
evaluated at 5. The surface stress tensor is also assumed 
to be “Newtonian” and defined by a Boussinesq-Scriven 
constitutive relationship of the form [4, 12, 17] 

n* = al s + 2p $ E s + k J,(V* • V s ). (6) 


In this expression, rr, p s , and k s respectively refer to 
interfacial tension, surface shear viscosity, and surface 
dilatational viscosity. Also, l s , = I - riri is the surface 
unit tensor, V s = I s • V is the surface gradient, and 
E* is the symmetric and traceless surface rate-of-strain 
tensor defined by 

E s = h(V 5 v s ) ■ I, + I, • (V 4 -v' , ) t ] - \l s (V s ■ v s )- 

(7) 

The derivation of the mechanical energy equation 
begins by dot multiplying the momentum equation in 
(1) by v, the momentum equation in (2) by v, integrat- 
ing over the respective material volumes Y{t) and I "(f), 
and adding the resulting equations. Assuming that the 
concentration of surfactants T varies only slightly from 
the equilibrium concentration T a , the surface tension 
may be modeled as a linearly decreasing function of the 
form 

a(T) =a 0 + /?( r - r 0 ) , /3= ^(T 0 ) < 0 . (8) 

With the aid of the boundary conditions (5), as well 
as the bulk and surface Reynolds Transport Theorems 
and Divergence Theorems [2, 12], the following energy 
equation is finally obtained [ 1 3] 


f t {K.B. + S.E.} = 


/ 2//(E : E) dV - [ 2/7 (E : E)dV 

J i (0 J\'(t) 


/ [2/i s (E 4 : E s ) + k 5 (V 4 • v*) 2 ] dS 

Js(t) 


-[ j3{T - r o )(V s ■ V s ) dS, 
Js{t) 


(9) 
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3 BOUNDARY INTEGRAL METHOD 


where 


be obtained to describe the time- dependent dynamics: 


K.E. + S.E. = 



+ / a 0 dS (10) 

Js(t) 

is the total kinetic energy plus surface potential energy 
of the system. The first two terms on the right-hand 
side of (9) represent dissipation in the two bulk fluids. 
The next two terms are similarly identified as dissipa- 
tion terms arising from surface shear and dilatational 
viscosities. The last term on the right-hand side of (9) 
requires further attention. This term can be shown [13] 
to be either dissipative or provide an additional surface 
energy storage depending on the surface Peclet number 
characterizing surfactant transport. The complete sur- 
factant transport equation for an insoluble surfactant is 
itself given by 

% + v s ■ v s r + (v s • v s )r = v s • (D s v a v), an 

at 

where D s is the surface diffusivity of surfactants. 

The mechanical energy equation obtained above may 
be used to derive an approximate ordinary differential 
equation (ODE) which describes the evolution of pure 
shape modes of three-dimensional axisymmetric liquid 
drops [13]. The key results are outlined here and the 
reader is referred to the original article [13] for details. 
The approximation involves using the potential flow so- 
lution to estimate the kinetic energy and the dissipation 
integrals in the bulk, the oscillatory Stokes boundary 
layer velocity field to estimate the dissipation rate in 
the thin layers surrounding the interface, and the sur- 
face velocity from this analysis to estimate the surface 
dissipation integrals. For the quadrupole mode of oscil- 
lation, with the radial coordinate of surface points given 
by 

r = a o [l + f(t)P 2 (cos0)], (12) 

and the surfactant concentration by 

1(6, t) = I o + g(t)P 2 (cos0), (13) 

where a Q and F 0 are the respective equilibrium radius 
and surfactant concentration, f(t) and g(t) the time- 
dependent amplitudes of deformation and concentration 
change (both assumed small), P 2 the second Legendre 
polynomial, and 6 the polar angle measured from the 
axis of symmetry, the following system of ODEs may 




f + 4 o 0 f = -$g 


5 fia 0 + I2p s + k s + 


25ag 

12\/2 


y/W2,oHP 




m + ^9(t) = Tom. 

”0 

Here, u 2 , 0 is the base frequency of the quadrupole shape 
mode given by (8cr 0 /pa^) I/2 . The coefficient of / on 
the right-hand side of the first equation represents damp- 
ing due to viscous dissipation in the bulk liquid inside 
the drop, due to surface shear and dilatational viscosi- 
ties, and due to the Stokes boundary layer in the gas 
surrounding the drop. In addition, the term which cou- 
ples that equation to surfactant concentration g(t) can 
be partially dissipative, depending on the value of the 
surface Peclet number D s [13]. These results 

are obtained in the limit where the density and viscosity 
of the drop phase are large compared with those of the 
surrounding fluid. 


3 Boundary Integral Method 

The numerical problem considers the shape oscillations 
of a /wodimensional inviscid drop in vacuum without 
gravity. It is assumed that the drop has an undeformcd 
equilibrium radius a 0 and the interface is highly con- 
taminated with insoluble surfactants so that the interfa- 
cial properties & 0 , // s , and k s are constants. Parameter- 
izing all the variables with arclength s, using the invis- 
cid stress tensors II = —pi and IT = 0, and expressing 
the gradient V = s d/ds + nd/du and velocity of the 
interface v s = si' s + hv u in terms of local coordinates 
s and ri, the tangential and normal components of the 
dynamic boundary condition in (5) take the respective 
forms: 

° = <,4) 

V = C{„„ + + »„(?]}, (15) 

where C is the local curvature of the two-dimensional in- 
terface. By first noting that v 9 is necessarily periodic in 
total arclength L and integrating (14) twice with respect 
to arclength around the drop these two components may 
be combined to obtain 

p = C{<r 0 + K s B{t)}, (16) 
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4 CONCLUSIONS 


where the time-dependent quantity B(t ) is defined by 


£?(f) — ^ J ds. 

(17) 

Nondimensionalizing time with (pi? 3 j a 0 ) ! , length 

with a 0 , and mass with pa 3 the boundary integral for- 
mulation of the governing equations for this system at a 
particular instant in time are 

,L 

/ k\(si,s)[q(s) - q{si)]ds = <£(«;), 

Jo 

(18) 

if>(si) = - f K- 2 {si,s)[q(s ) - q{si)]dsk, 

Jo 

09) 

!£(*)■ 5 |V(«)I’-C(«)U + «:*). 

(20) 


(21) 


Equation (18) represents the regularized double-layer 
potential boundary integral formulation of Laplace’s eq- 
uation for potential flow [141. I\\ is the weakly singular 
kernel defined as the projection of the gradient of the 
two-dimensional Green's function for Laplace's equa- 
tion in the direction normal to the interface, q(s) is a 
distribution of dipole densities around the drop inter- 
face, and (j> is the scalar velocity potential, tp = t k is 
the vector velocity potential perpendicular to the plane 
and is related to the velocity of the surface through its 
curl, Vx^=v*. This vector velocity potential also re- 
lates to the distribution of dipole densities through an in- 
tegral containing the weakly singular kernel /T>, which 
is defined as the cross product of the unit normal with 
the gradient of the two-dimensional Green's function. 

The dynamic boundary condition (16) has been used 
with Bernoulli's unsteady equation for pressure to ob- 
tain (20), where the time derivative is with respect to an 
observer moving with the velocity of the interface v s 
and k* - k 8 {<t 0 / pi? 3 ) 1/2 is the nondimensional sur- 
face dilatalional viscosity. Equation (21) represents the 
kinematic boundary condition. 

The above equations were discretized by dividing 
the periodic boundary into N equally-spaced nodes in 
the interval 0 < s < L. All derivatives were calculated 
using standard @[( A.s') fi ] finite-difference schemes and 
the regularized integral relations in (18) and (19) were 
discretized into matrix relations using a trapezoidal qu- 
adrature rule between the nodes. The numerical scheme 
first initialized the shape of the drop x and the scalar ve- 
locity potential <j>. The velocity of the interface was then 
calculated using LU decomposition to solve the matrix 


equivalent of (18), performing the matrix multiplication 
in (19), and taking derivatives of the vector velocity po- 
tential with respect to arclength. Using the updated ve- 
locity of the interface and curvature, the scalar veloc- 
ity potential and drop shape could be integrated in time 
using a fourth-order Runge-Kutta scheme. To prevent 
clustering and allow for the calculation of the deriva- 
tives, the nodes were redistributed to equal spacing in 
arclength after each time-step. The accuracy of the nu- 
merical method was checked by calculating the con- 
served quantities of total energy and volume in time. 
Interestingly, the numerical problems reported by [91, 
arising from the instability of modes with wavelength 
twice the nodal spacing, did not appear in these two- 
dimensional calculations. 

Figure 1 shows an example of the damping effects 
of surface dilatational viscosity for an initially perturbed 
two-dimensional inviscid drop. The calculation used 40 
nodes to simulate nearly 40 oscillation periods in the 
quadrupole mode of moderate initial amplitude. Al- 
ternate plots of similar energy versus time curves re- 
veal that the attenuation in time due to surface dilata- 
tional viscosity in two-dimensional drops cannot be rep- 
resented by an exponential or power law. 


4 Conclusions 

An energy equation has been derived for the general 
case of a viscous drop suspended in a viscous medium 
with surfactants contaminating the interface. It contains 
terms clearly identifying dissipation contributions from 
the viscous effects in the bulk fluids, surface shear and 
dilatational viscosity effects in the interface, and surfac- 
tant transport. 

An efficient numerical boundary integral method has 
been developed which incorporates the effects of a con- 
stant surface dilatational viscosity in simulations of an 
oscillating two-dimensional inviscid drop. Surface di- 
latational viscosity is shown to have a significant damp- 
ing effect on the otherwise undamped inviscid oscilla- 
tions. This damping was found to be neither an expo- 
nential nor a power law. 
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Figure 1: Total energy of the drop versus time after an initial condition r - 1 + 0.1 cos(20) and ^ - 0 for the cases 
= 0 (solid line), and k* = 1 (dashed). 
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